Calibration of contact parameters of sandy soil for planting tiger nut based on non-linear tools

A methodology combining physical experiments with simulation was employed to acquire contact parameters of sandy soil precisely for planting tiger nuts in the desert area of Xinjiang. The stacking angle under different parameter combinations was applied as a response value. Through the Plackett–Burman test, several factors that have a significant influence were determined. The steepest ascent test was conducted to establish the finest scope of values for these parameters. The stacking angle was considered the response variable, and non-linear tools were used to optimize these parameters for simulation. The findings showed that applying response surface methodology (RSM) resulted in a relative error of 1.24%. In the case of BP-GA, the relative error compared to the physical test value was 0.34%, while for BP, it was 2.18%. After optimization using Wavelet Neural Network (WNN), the relative error was reduced to only 0.15%. Results suggest that WNN outperforms the RSM model, and the sandy soil model and parameters generated using WNN can be effectively utilized for discrete element simulation research.

www.nature.com/scientificreports/deposition angle was determined to be 26.586°,exhibiting a relative error of 1.22%, thus confirming the accuracy of the simulation outcomes.Wang 7 investigated the calibration of contact parameters among sand particles, with findings indicating that the variance between virtual simulation and actual testing was less than 5.1%.Wu 8 , Lin 9 , and Wang 10 analyzed and calibrated the contact parameters of manure.All three studies confirmed the viability of the methodologies, thereby yielding commendable experimental outcomes.Coetzee 11 and Pasha 12 demonstrated that particle models of various shapes necessitated parameters during the simulation, which should be calibrated via experimentation.Moreover, Liu, Wang, and Liu [13][14][15] calibrated contact parameters of crops such as wheat, corn, and mini-potatoes, respectively, providing valuable references for parameter optimization and performance analysis of crop harvesting mechanization.The above studies have all adopted the discrete element calibration method, which allows for relative motion between individual cells, does not necessarily meet the conditions of continuous displacement and deformation coordination, has fast calculation speed, and requires small storage space to solve problems with large displacement and nonlinearity, and has achieved promising experimental results.
Furthermore, machine learning techniques have found extensive applications in engineering.For instance, references [16][17][18] demonstrate modeling using the cosine amplitude method, and relevant algorithms were then employed to analyze and predict various engineering issues.Simultaneously, numerous machine learning studies to solve engineering problems have proposed several valuable and innovative designs [19][20][21][22] .
The studies mentioned above have yielded promising results; however, there is still scope for further improvement.Based on the study mentioned above, this research focuses on sandy soil found in the desert areas of Xinjiang, where tiger nuts are cultivated.A combined approach was utilized to calibrate the contact parameters of sandy soil's discrete elements 23 .Afterward, the obtained crucial parameters were optimized and fitted by applying non-linear tools, followed by the acquisition and comparison of predicted values.Consequently, a model was developed, serving as a conceptual background for the design of devices in the mechanized production of economic crops like tiger nuts.

Materials
The sandy soil was acquired from a tiger nut plantation in the 54th Regiment of the Third Division in Xinjiang.The 54th Regiment is in Shache County, on the southwestern frontier of Xinjiang, at the northern foothills of the Kunlun Mountains, and south of the Pamir Plateau.It lies within the mid-upper reaches of the Yarkand River alluvial fan plain, positioned between the Taklamakan Desert and the Bugur Desert, with the terrain sloping from west to east.The geographical coordinates are approximately between 76°1′57″ to 77°46′30″ E longitude and 37°27′30″ to 39°15ʹ N latitude.A bulk density meter was employed to determine its density.The standard for the determination of the particle size distribution is GB/T 50123-2019.The particle size distribution was determined by screening soil by applying a standard particle size sieve.Additionally, moisture content was measured, as presented in Table 1.The intrinsic parameters for sandy soil and steel plate were determined based on previous studies 6, [24][25][26] and the GEMM database of EDEM, as listed in Table 2.
As an intrinsic characteristic of bulk materials, stacking angle can provide insights into their flow and friction properties [27][28][29] .A combinatorial approach 11 was employed to calibrate parameters.Injection methodology was utilized to measure the actual stacking angle.Plackett-Burman design, conducted using Design-Expert 8.06, was applied for a multi-factor significance screening test and analysis to identify parameters that have a crucial influence on stacking angle.Box-Behnken design RSM was applied to form and ameliorate a regression model.Compared to the central composite design, the Box-Behnken design requires fewer trials and typically features a more uniform layout of test points.Additionally, it excludes pivot points on the boundaries of the test range, thereby reducing the extreme conditions of the experiment.Therefore, this study opted for the Box-Behnken design.The actual stacking angle was adopted as the objective value.A regression model was also used to determine the combination of crucial parameters.Subsequently, a simulation was performed by applying optimal parameters.The simulated value was contrasted with the actual value to validate the validity of the calibrated parameters.

Contact model betwixt particles
The hertz-Mindlin (no slip) contact model is known for its accuracy in force calculation.This model proves highly efficient given the minimal cohesive action betwixt particles in sandy soil 33 .
Assuming that two particles with radii R 1 and R 2 are brought into elastic contact, the mathematical calculation for normal overlap α is r 1 and r 2 are sphere center position vector.Contact radius of contact surface is a.R * is equivalent particle radius, which is computed as follow, Normal force betwixt particles F n is E * is equivalent elastic modulus, which is computed as follow, (1) υ rel τ means tangential component of relative velocity.Rolling friction be described as, where, µ r means rolling friction coefficient, R i represents distance from contact point to center of mass.In addition, ω i is unit angular vector of object at contact point.

Establishment of geometric model
To establish a suitable particle model, the mass distribution of sandy soil particles ranging from 0 to 0.25 mm was scrutinized, and it was found that they accounted for 92.29% of the total mass among the selected test samples, making them representative of the overall sand sample.Therefore, these particles were selected for modeling purposes.It should be noted that the size and shape of the sand particle impact the results and the computational efficiency of the computer.Considering laboratory conditions and physical properties of sand particles, spherical particles with a size of 1 mm were selected.
Using a physical test model as a basis, a model was created in SolidWorks software, which was integrated into EDEM as a simulation model for the funnel device, utilizing the Hertz-Mindlin (no slip) contact model.Figures 3 and 4 display particle model and funnel model, respectively. (5) It is determined that the initial radius of particles generated in the simulation procedure should be 1 mm, and they were randomly generated.To prevent the generation of tiny particles, the radius of the generated particles was limited to between 0.5 and 1.25 times the initial radius.
Dynamic generation of particles in the funnel was performed using a virtual particle factory above the funnel.The overall weight produced was 0.45 kg, with a rate of 0.2 kg/s.The time interval for data storage was established as 0.01 s, while the fixed time step was 20% of the Rayleigh time step.Furthermore, the mesh size was three times the smallest spherical element size 6,34,35 .Values of the other parameters are shown in Table 3.

RSM test methodology
Plackett Burman design (PBD) screening experiment.A PBD screening experiment was adopted to ascertain the importance of every individual factor.The test examines the correlation between objective response and every individual factor, contrasting the difference between the two levels of every individual factor.The stacking angle of sandy soil was selected as the response value for screening contact parameters.Two levels were chosen for every individual of the six factors.Simulation parameters were determined based on characteristics of sandy soil used for planting economic crops in desert areas of Xinjiang.Table 4 shows levels of factors, along with corresponding symbols.
Steepest ascent test.The steepest ascent test is a methodology that can quickly identify the area with optimal response value while minimizing the number of tests required.It was conducted by increasing or decreasing values of significant parameters according to predetermined step length 36 .For residual parameters, the median  where, σ is stacking angle of simulation test (°).θ is stacking angle of physical test (°).
Box Behnken design response surface experiment.Based on the principle of response surface design 37 , the Box-Behnken test was conducted with three levels of significant parameters (low, medium, and high).
Additionally, five central points were included to estimate the error.Factor values are provided in Table 5. Factors corresponding to each symbol are consistent with Table 4.

Machine learning algorithm
The Back Propagation (BP) algorithm is a feedforward learning algorithm that utilizes a multi-layer perceptron, also known as BP neural network, to perform backpropagation.This algorithm is iterative and involves forward propagation, backpropagation, and parameter updates.On the other hand, Genetic Algorithms (GA) offer advantages such as efficiency, parallelism, and global search.Normalization processing must be performed before using the BP-GA integrated optimization algorithm for data training.This is followed by selecting one feature value as the output and the remaining feature values as the input.Upon completion of training, the BP-GA integrated optimization algorithm generates prediction results.Wavelet Neural Network (WNN) is constructed based on the wavelet transform theory, and its principle is similar to that of the Backpropagation Neural Network (BPNN).The primary characteristic of WNN is that its hidden layer neuron activation function is a wavelet basis function, which effectively utilizes the localization property of wavelet transform and the large-scale data parallel processing and self-learning capability of neural networks.Consequently, it exhibits a vital ability for approximation and achieves fast convergence speed.The specific flowchart is shown in Fig. 5.
Train and optimize stacking angle using the WNN optimization approach.Before commencing the experiment, three sets of data were randomly chosen from every amalgamation of factor levels, leading to a total of 51 data sets being extracted.WNN was employed to learn and train these 51 data sets as training sets.Subsequently, the algorithm's performance in predicting stacking angle was assessed using 17 other test sets, which can be seen in the table.
Applying WNN to predict the stacking angle response value involves predicting the dependent variable by considering the explanatory variables.Previous experiments have identified B, D, and E as the three most significant influencing factors, and this experiment employs them as explanatory variables.X 1 , X 2 , and X 3 were used as input, and Y was used as the output (as shown in the topology diagram in Fig. 6).B, D, and E were represented by the numbers 1, 2, and 3, respectively.To enhance training outcomes, it is essential to normalize the training sample set in the preprocessing stage, considering the significant disparities in the magnitude of each input variable.The objective of this normalization process is to attain superior training results.

Simulation experiment and result analysis PBD screening test
Table 6 displays the design and outcomes of PBD.ANOVA analysis was conducted on outcomes applying Design Expert 8.06 38 , and the significance of parameters is displayed in Table 5, which reveals that factors A, B, and D exhibit a positive repercussion, indicating that stacking angle enhances as these factors enhance.Conversely, factors C, E, and F have a negative effect.The significance sequence of every factor was ascertained based on the P value of factors shown in Table 7.Three factors, E, D, and B, which have a higher significance, were selected in the steepest ascent experiment and the Box-Behnken design experiment.

Steepest climb test result analysis
Findings from the steepest ascent test are displayed in Table 8.Findings revealed that the stacking angle gradually enhanced as factors B, D, and E lessened.Furthermore, the relative error between the actual stacking and initial angles was lessened initially and then enhanced.The most minor relative error was observed in the No. 5

RSM analysis test
The experimental plan comprises 17 test points, encompassing 12 analysis elements and 5 zero estimation errors.9, indicate that three significant parameters investigated, i.e., B, D, and E. Design Expert 8.06, were applied to find a quadratic regression equation.The equation is as follows, ANOVA results are displayed in Table 10.For fitting the mathematical model, P = 0.0003, revealing an excellent fitting degree.P value of B was less than 0.01, P values of D, E, interaction item B × D, and interaction item B × E were all less than 0.05, confirming their significance on stacking angle and validating regression.For Lack of Fit, P = 0.1529 > 0.05, revealing a tiny proportion of atypical errors in fitting betwixt obtained regression equation and actual value, with a splendid fitting degree without lack-of-fit phenomenon.The coefficient of variation CV = 1.7% indicated high validity.The coefficient of determination R 2 = 0.964 and adjusted coefficient of determination R 2adj = 0.918, all-surpassing 0.9, demonstrating that the proposed model precisely represents an actual situation.Additionally, Adep Precision = 14.948, demonstrating the proposed model's high precision.
According to the results in Table 11, factors that have no significant impact (D × E, B 2 , D 2 , and E 2 ) were eliminated.ANOVA results of the modified mathematical model are displayed in Table 11.For fitting mathematical model, P < 0.001; for Lack of Fit, P = 0.1043; coefficient of variation CV = 2.08%; determination coefficient R 2 = 0.915; adjusted determination coefficient R 2adj = 0.877; test precision Adep Precision = 15.969.It reveals that the proposed model had a splendid fit, validity, and accuracy and had specific improvements that were in contrast with those before optimization.After optimization, the regression formula is,

Analysis of interaction effects
Based on the ANOVA results of the regression model, interaction terms involving B and D and B and E had a noteworthy repercussion on stacking angle.The response surface slope's steepness further confirms the significance of these factor interactions.Moreover, the contour lines' shape indicated the strength of this interaction 39 .Response surface outcomes for the stacking angle of sandy soil are illustrated in Fig. 7.
When set E at 0.3, the response surface plot for the interaction between B and D is presented in Fig. 7a,b.It could be observed that the response surface curve for B was steeper contrasted to the increasing direction of D. The Contour density of B was higher than the density along the increasing direction of D, revealing that B had a crucial repercussion on the stacking angle.The contour shape appeared to be ellipsoidal, suggesting a notable interaction betwixt B and D. When D was held at a certain value, the stacking angle enhanced as B enhanced.Similarly, when set B at a specific value, the angle tended to enhance as D enhanced.Notably, when B took on an immense value, the angle showed a significant enhancement as D enhanced.
When set D at 0.5, the response surface plot for interaction between B and E is presented in Fig. 7c,d.The response surface curve for B was steeper than E's increasing direction.The Contour density of B was higher than the density along the increasing direction of E, revealing that B had a crucial repercussion on the stacking angle.The contour shape appeared to be ellipsoidal, suggesting a notable interaction betwixt B and E. When E was held at a specific value, the stacking angle enhanced as B enhanced.Similarly, when set B at a specific value, the angle tended to lessen as E enhanced.Notably, when B took on a significant value, the angle significantly lessened as E enhanced.

Comparison of machine learning models
This paper suggests conducting a comparative analysis of different machine learning algorithms to determine the most efficient fitting algorithm.Table 10 presents a comparative analysis of the non-linear algorithms.Figure 8a presents a visualization of WNN training residuals.The results demonstrate that the WNN surpasses the others in achieving lower relative errors.Besides, the BP algorithm registered a mean absolute error (MAE) of 0.28 and a root mean square error (RMSE) 0.64.Meanwhile, the BP-GA algorithm showed improved performance with an MAE of 0.11 and an RMSE of 0.18.WNN algorithm outshone both, with the lowest recorded MAE of 0.068 and an RMSE of 0.035.Figure 8b illustrates the results of the train set, while Fig. 8c,d depict the RMSE of the training process and MAE of the training process, respectively.This outcome underscores the substantial limitations inherent in both BP and BP-GA algorithms about their precision in identifying optimal weights and thresholds.This excessive workload adversely impacts the accuracy and stability of its predictive outcomes.
This experiment has identified WNN as the most excellent optimization machine learning algorithm.Through analysis, it becomes evident that the training residual gradually decreases as the epoch progresses.However, the RMSE and MAE exhibit significant fluctuations throughout the training process.

Comparison of non-linear tools
RSM parameter optimization Design Expert was harnessed to determine the finest values of three significant factors by solving an optimized regression equation.The optimal value of B was 0.63, D was 0.55, and E was 0.43.Median values of nonsignificant parameters were chosen (A was 0.4, C was 0.2, and F was 0.22).The abovementioned parameters were adopted to validate the accuracy of optimal parameter amalgamation, and other parameter configurations remained unchanged.EDEM was applied to execute simulation tests on stacking angle.The average value acquired from three repeated simulation tests was 31.75°.Its relative error from the actual value of 32.15° was 1.24%.Optimization constraints for solving it are defined as,

Figure 8 .
Figure 8.The results of train set.

Table 1 .
Fundamental elements of sandy soil.

Table 2 .
Intrinsic parameters of sandy soil and steel plate.

s ratio Shear modulus Density(g/cm −3 )
www.nature.com/scientificreports/where,E 1 represents the elastic modulus of Particle 1, ν 1 denotes the Poisson's ratio of Particle 1, E 2 represents the elastic modulus of Particle 2, and ν 2 denotes the Poisson's ratio ofParticle 2.Normal damping force can be acquired by later calculation.

Table 4 .
PBD. values were selected.Specifically,A was 0.4, C was 0.2, and F was 0.22.The test scheme can be found in Table7, and the relative error (N, %) is shown in the formula(14).

Table 6 .
Scheme and results of PBD.G, H, J, K and L indicates blank column.

Table 8 .
Scheme and results of the steepest ascent experiment.

Table 9 .
Design and results of Box-Behnken.

Table 11 .
ANOVA of modified model of Box-Behnken experiment.